#### Age-based yield data from Arak (1967) and checks

library(ggplot2)

## Configurations

do.figures <- F

years <- 100
mciters <- 100000

## Yield data from Arak (1967)

yage.age <- c(1:3, 5, 8, 11, 14, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 55.5, 65.5, 75.5, 100)
yage.bar <- c(0, 0, 99, 441, 634, 541, 536, 499, 392, 435, 405, 357, 365, 355, 406, 440, 361, 299)

yage.spanlo <- c(3, 4, 7, 10, 13, 16, 21, 26, 31, 36, 41, 46, 51, 61, 71, 81) - 1
yage.spanhi <- c(3, 6, 9, 12, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 100)

yage.weight <- c(rep(1000, 3), rep(100/3, 4), rep(10/5, 7), rep(1/10, 3), 1/20)
yage.int <- predict(smooth.spline(yage.age, yage.bar, w=yage.weight, df=length(yage.age) - 5), 1:100)$y
yage.int[1:2] = 0

ybar <- mean(yage.int)

if (do.figures) {
    ## Yield data
    ggplot(data.frame(year=1:100, yield=yage.int), aes(x=year, y=yield)) +
        geom_line() +
        geom_errorbarh(data=data.frame(yage.age=yage.age[3:length(yage.age)], yage.spanlo, yage.spanhi, yage.bar=yage.bar[3:length(yage.age)]),
                       aes(x=yage.age, xmin=yage.spanlo, xmax=yage.spanhi, y=yage.bar), colour=2) +
        theme_bw() + ylab("Yield (kg / Ha)") + xlab("Plant age (years)") + scale_x_continuous(expand=c(.01, 0))

    ## Repeated shape under management
    ggplot(data.frame(year=1:100, yield=c(yage.int[1:33], yage.int[1:33], yage.int[1:34]), yield2=c(rep(0, 4), yage.int[5:33], rep(0, 4), yage.int[5:33], rep(0, 4), yage.int[5:34])), aes(x=year, y=yield)) +
        geom_line(linetype=2) + geom_line(aes(y=yield2)) + geom_hline(yintercept=395, col=2) +
        theme_bw() + ylab("Yield (kg / Ha)") + xlab("Plant age (years)") + scale_x_continuous(expand=c(.01, 0))
}
